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1. Introduction 



Modeling time-series with linear pole-zero AutoRegressive-Moving Average (ARMA) 
models has numerous applications in signal processing. This problem is in general non linear and 
most ARMA modeling techniques are iterative in nature. The Iterative Prefiltering (IP) method 
has the advantage of computing potential non-minimum phase representations which may be useful 
in time-domain modeling. The original IP minimization procedure is an ill-conditioned problem 
which has classically been solved using a least-squares approach. This work presents a modification 
of the classical IP technique in which the least-squares iteration step is replaced by a Total Least 
Squares (TLS) step to take advantage of the statistical properties of the TLS method. Results show 
that improvements in the modeling performances may be obtained by using the TLS-based IP 
method when modeling signals distorted by white Gaussian noise. 

Section 2 presents a review of the TLS method, and illustrates its main advantages as 
compared to the LS method. Next, Section 3 presents the classical Iterative PrefiJtering method and 
the proposed TLS-based IP method. Finally, conclusions are given in Section 4. 

2. Total Least Squares Problem 

2.1 Introduction 

The Total Least Squares method (TLS) has been introduced recently as a substitute to the 
classical Least Squares (LS) method for solving overdetermined linear systems of equations when 
all data involved in the computation are corrupted with noise (errors) [1,2]. Consider the problem 
of solving the overdetermined system of equations Ax = b, where A has dimension m*n and b has 
dimension m*l. The LS solution x, s is obtained by minimizing the Euclidean norm ||Ax-b|| 2 . 
Solving for the LS solution is equivalent to solving 
the following minimization problem: 



minimize \\h~b!. || 2 
subject to bL e R(A) 



The above equation is satisfied when b’ is the orthogonal projection of b onto R(A). Thus, the LS 
vector b’ can be viewed as the perturbation of b so that b’ can be generated by the range (or the 
column space) of A. As a consequence, the LS solution assumes that errors can occur only in the 
vector b and not in the matrix A. However, this assumption is not always realistic as errors such 
as sampling, modeling, and instrumentation errors may force inaccuracies in A also. The TLS takes 
into account the fact that errors may occur both in the matrix A and in the vector b and solves the 
following minimization problem: 

minimize^ £ e W{A\h\-[A\b ] || 2 

subject to b e R(A) 

Thus, the TLS vector B can be viewed as resulting from the perturbation both from the columns of 
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A and the vector b so that 6 belongs to R(A). An illustration of the differences between the LS and 
the TLS solution for a system with column space of dimension two is shown in Figure 1. 

22 SVD-Based Solution of the TLS Problem 



The Singular Value Decomposition (SVD) is used to solve the TLS problem [1,2]. Recall 
that the system to be solved is of the form Ax = b, which can be reformulated as: 



m 



-i 



=q. 



The TLS solution is found from the SVD of the augmented matrix 

C = [A\b] = 1/2 K w . 



The solution is said to be generic or non generic depending on the numerical characteristics of the 
eigenvalue matrix E and the right singular vector matrix V. 

22.a TLS Generic Solution 

The TLS solution x-tls 1S unique and said to be generic if the singular value o n >o n + 1 and 
V n+i.n+i*0. In such a case the solution is given by: 

Xtls = • 1 / v n + i.n + i[ v i,„ + i-A’ n ,„ + i] 1 - 

When the p smallest singular values are equal, i.e., o n . p >o n . p+1 = ... = o n+1 , the TLS solution 
to the linear system is not unique, and any linear combination of the singular vectors associated with 
the multiple minimum singular value can be chosen provided that it is normalized properly. The 
resulting minimum norm solution can be shown to be equal to [1,3] 

V 2 Q V { g 



Matrices Vj and V 2 are submatrices of the matrices of the right singular vector matrices V; and V 2 , 
where V 2 is the matrix of right singular vectors associated with the multiple minimum singular value, 
and Vj contained all the other right singular vectors. V, and V 2 are defined as: 



^2 = 



c' 


, and V. = 


g l 


A 


5 1 


A 



where the vectors c‘ and g 1 respectively represent the first row of V 2 and Vj. 

Again, the TLS solution exists only if all v nt]1 *0 for i = n-p+ l,...,n+ 1. Van Huffel showed 
that no generic solution exists only when the matrix A is nearly rank deficient, or when the set of 
equations is conflicting [4], The set of equations is said to be conflicting when o n . p+I =...=o n + 1 are 
large. In such a case, trying to model the data using a linear model is inaccurate, and a better 
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option should be chosen by the user. When o n . p+I =...=o n+1 are small, two options are possible. The 
first option for the user is to remove dependent columns of A to reset the problem so that the 
resulting modified A matrix is regular, and to solve a generic TLS problem. How to pick such 
columns is addressed in [2,14], The second option is to solve a nongeneric TLS problem which is 
obtained from adding additional constrains on the generic TLS problem in order for it to be 
solvable. 

22.b The Nongeneric TLS Problem 

The nongeneric TLS problem is addressed in details in Van Huffel et al. [4,5,7], The 
resulting solution is given by: 



when o„. p > o n . p+1 



—TLS = -V V „ + l.„.p[ V l.„.p>--->V„,„. p ]‘, 

° n+ i> v „ + u = ° for i= n ’P+ 1 -- n+ 1, and v n+ln . p # 0. 



23 Applications of the TLS to Signal Processing Problems 

Numerous researchers in Signal Processing have reformulated problems in terms of the TLS 
technique. Applications can be found in array processing [11], in system modeling [10], and in 
frequency estimation of sinusoidal signals [9,10-13] for example. Results show that improvements 
in the performances can be obtained when reformulating LS problems in a TLS set-up for the 
applications mentioned above. This result motivated our investigation of the TLS technique to 
time-series data modeling. ARMA modeling of time-series is a non-linear problem which has been 
extensively studied in the literature [17], The Iterative Prefiltering (IP) method is an iterative 
linearized formulation which was originally proposed by Steiglitz and McBride to identify linear 
system transfer functions [15,16], and the IP method can be reformulated to model time-series data 
distorted with noise. Note that the method is not insured to converge to the optimal solution and 
De Moor showed that it converges to a suboptimal solution only [18]. Nevertheless, simulation 
results have shown that useful models can still be obtained using the IP method [19]. 

3. The Iterative Prefiltering Method 

3.1 Introduction 



The Iterative Prefiltering method attempts to model a time-series data with an ARMA(P,Q) 
system using a time-domain approach by minimizing the error between the data and the impulse 
response of the system to be estimated. Using the Z-domain, the model function is given by: 



H(z) = 



A(z)' 



where 
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( 1 ) 



p p 

>4(z) = Y, a(k)z~ k and B(z ) = Y b(k)z ~ k . 
jU o k = o 



The error function for the problem is defined by 



E(z) = X(z)~ 



B(z) 

>4(z) 



X(zM(z)-l*fc) 

>4(2) 



Non-linearity in the problem is due to the denominator A(z) in the error function E(z). The 
problem may be linearized by replacing the error function given in eq. (1) with the iterative error 
function 



£(0(z) = X(zM (0 fc)-B <0 (z) 

A ^ _l) (z) 



where indices (i) and (i-1) refer to iterations (i) and (i-1). At each iteration the quantities A 0) (z) 
and B (l) (z) are chosen to minimize the error E (0 (z). The estimation of the ARMA transfer function 
is done in the time domain using the time series x(n), for n = 0,...,N-l, by rewriting E (1) (z) in the 
time-domain as 

E«\ n ) =X <W(„) *gi i >-h m (n\b i0 . n= ( 2 ) 



v-'ith: 



£ w (n) n=0, ,N - 1 



h^(n)=[h°\n)„ ,h°\n-Q)) T 

a^=[l,a ( ‘\l), y\P)] 
blH=[b w (0),...,b l, \Q)) 



where h <0 (n) is the impulse response of H <0 (z)= 1/A (ll (z), and x (1> (n) represents the output of x(n) 
through the filter H 0) (z). Eq. (2) expressed in matrix form becomes: 





0 


;c ( °(0) 




a u \ P) 




/i w (0) 


0 


0 




b ( ‘\ 0) 


= 


0 


x (0( 0 ) *<0(1) 




a ( '\P- 1) 


- 


A ro (l) 


h u \0) . 


0 




b (, \ 1) 




x u \N-P) . 


.. x u \N-2) x (, \N-l) 




1 




h°\N- 1) 




.. h°\N-Q) 




b°\Q) 



which may be rewritten as: 
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E (0 = X (i) dH-H li) til. 



(4) 



Equation (4) may be decomposed into two parts using the orthonormal projection matrix P onto 
the column space of X (l) , as noted by McClellan and Lee [17], which leads to: 

El 2 = 2 ) 

(5) 

with P = 



Recall that the IP iteration minimizes the expression ||E (1) || 2 , which leads to the following 
minimization problem [17]: 

• (a) minimize p-X^a® subject to a w (0)=l, 

• (b) solve 

Step (a) above may be rewritten as: 

[C,-d][a (1) (P),...,l] T = 0, 

where 

C = P i X (i) [l:N,l:P], d = P x X (i) [l:N,P+ 1], 

Step (a) of the IP method has classically solved using the classical LS approach, which 
accounts for errors in d only. However, errors occur both in C and d. A better fit of the data can 
be obtained by taking into account for errors in both in C and d, as the TLS set-up allows us to do. 
Thus, the TLS-Based IP iteration solves for (a) using a TLS approach [20], which leads to the 
following minimization problem: 

• (a) Minimize P x X (,) a (1> subject to a (,) (0) = 1, using a TLS method. 

• (b) solve ^=[H (W H {i) y x H (i)lt X^. 

3 2 Data Scaling 
3.2.a Introduction 

Van Huffel has shown that when errors in the TLS matrix and right-hand-side vector are 
uncorrelated with zero-mean and equal variance, then under mild conditions the TLS solution is 
a strongly consistent estimate of the true solution of the unperturbed system [4]. For the problem 
considered, errors in the data result from errors in the signal x(n) to be modeled. Assuming the 
data to be modeled is distorted with Gaussian white noise, the errors in H (l) are correlated as they 
are obtained as the output of an AR linear system. However, consistency in the TLS estimate may 
still be retained by using the Generalized TLS (GTLS), as proposed originally by Van Huffel 
[2,4,17], Thus, data scaling is needed to insure diagonal error covariance matrices, and the problem 
becomes to estimate the error covariance matrices F (,) and G (l) at each iteration given by: 



F (i) = E[N H(i) N (i) ], 
G <i) = E[N (i) N (i)H ], 
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where N (,) represents the errors (noise contribution) in [C,-d] at iteration (i). 

3.2.b Computation of the Matrix F"’ 

Assume that the noisy signal x(n) is defined by x(n) = s(n) + w(n), where w(n) is wss white 
noise, and s(n) is an ARMA signal. The noise correlation matrix F 0) may be expressed in terms 
of the correlation matrix R w (1) obtained from the AR process 1/A (l| (z): 

F (i) = E[N m N ( '■)] = E^P^'YiP^ 0 )] 

= E[W iW P jJ1 P L I0°] 

= E[W i ‘ )ll P x W i ' ) ] = (6) 

R<0 _ £« 



Computations show that the components of E w (l) in equation (6) can be expressed as: 

P q /Z(\Uq-j-l\), ( 7 ) 

<?A1 

where P t j is the (i,j) lh component of the projection matrix P defined earlier, and r w (l) (n) is the 
correlation function obtained when passing the white noise w(n) through the AR system with 
transfer function 1/A (1) (z). Using the fact that the noise distortion w(n) is assumed to be white, 
r w (l) ( 1 1 + q-j-h | ) is non zero only when l + q-j-k = 0. Thus, the double summation in equation (7) 
reduces to a single summation. Furthermore, the correlation function r w (0 (k) can easily be 
computed using the results by Dugre et. al. [6] who proposed an algorithm to generate a covariance 
sequence from its AR coefficients. 

3.2.C Computation of the matrix G"’ 



Again the error covariance matrix G (1) may be expressed in terms of the correlation matrix 
R w (l) obtained from the AR process 1/A l0 (z) obtained at iteration (i) by: 



F d) = E[N (0 N m ] = P^E[W i ' ) W iS)H ]]P x 
= P-R^P 1 



( 8 ) 



2>2.A GTLS Solution 

Applying the results presented by Van Huffel, the general TLS-based IP iterative solution 
for a (,) at iteration (i) is then given by solving: 



7 






“Q with RjP =chol(F (i) ), R ( £ = UZ'U H , 



with: 







rU) 




\ 


1 i 



E * =diag[Jo], ^0,...0] 



where r is the rank of R G (,) . Next, the estimate for b (l) can be obtained by replacing a (l) by its value 
in the expression: 



3.3 Simulation Results 

The TLS-based IP method is implemented on time-series data generated by an ARMA(3,4) 
in additive white noise and its performances compared with that of the LS-based IP method. In both 
cases, initial estimates for the iterative procedures are obtained using Yule Walker equations for 
the poles and Shank’s technique when solving for the zeros. The resulting noisy signal is modelled 
with an ARMA(5,6) using a sequence of length 60. The noise variance is chosen equal to 2 and the 
ARMA parameters are chosen equal to: 

a p = [1,-2,1.6725,-0.4613] 
b p =[l, -.9, 0.04, 0.3960, -0.29 12] 

Figure 2. a shows the noise-free signal, the noisy signal, and the estimated signal obtained 
during the first 5 iterations using the classic IP method. Power spectral estimates obtained for the 
noise-free signal, the noisy signal, and that obtained with the initial ARMA coefficients estimated 
using Yule Walker equations and Shank’s method, are shown in Figure 2.b. Figures 3.a and 3.b 
represent the same estimates obtained when using the proposed TLS-based IP scheme. Note that 
Figures 2 and 3 represent the best performances obtained during the first 5 iterations for both 
schemes. Figure 4 and 5 represent the true pole (a p ) and zero (b ) locations, and estimated poles 
(a ip ) and zeros (b jp ) locations obtained when using the IP and TLS-based IP schemes. This example 
illustrates the fact that the bias in the estimated pole locations obtained using the TLS-based 
technique is usually smaller than that obtained using the LS-based IP method when the SNR is 
medium to high. However, we noted that when the SNR is low, no distinct improvement is noted 
consistently when using the proposed technique. 

4. Conclusions 

This study has investigated the application of the TLS problem in the Iterative Prefiltering 
method for time-domain ARMA modeling. Results show that the pole bias is usually smaller when 
using the TLS-based scheme than when using the classical LS-based IP method, when SNR are 
medium to high. We note that, similarly to the LS-based method, the TLS-based IP method is not 
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TLS-based IP method. The main drawback in the proposed TLS-based IP method is the 
computational load increase. 
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Figure l.a) LS solution to Ax = b; b) TLS solution to Ax = b. 
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estimated signal: dash - noise free signal:solid - noisy signal; do t 




H_est (f ) ; dash - H_noise_f ree ( f ) : dot - H_itl (f) : solid 




Figure 2. a) Noise-free signal, noisy signal, LS-based IP estimated signal; b)Spectra: H_noise_free: 
noise-free signal, H_est:LS-based IP signal, H_itl; initial estimate. 



11 



estimated signal: dash - noise free signal;solid - noisy signal: do t 




H_est (f ) : dash - H_noise_f ree (f ) : dot - H_it 1 (f) : solid 




Figure 3. a) Noise-free signal, noisy signal, TLS-based IP estimated signal; b)Spectra: 
H noise_free: noise-free signal, H_est: TLS-based IP signal, H_itl: initial estimate. 
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a_ip; (* ) ap; (o) 




b_ip: (*) bp: (o) 




a_ip;(*) b_ip: (o) 




Figure 4.True pole (a p ) and zero (b p ) locations; LS-based IP estimated pole (a ip ) and zero (b ip ) 
locations. 
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a_ip: ( * ) ap: (o) a_ip;(*) b_ip: <o) 
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Figure 5.True pole (a p ) and zero (b p ) locations; TLS-based IP estimated pole (a ip ) and zero (b ip ) 
locations. 
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